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We study the rhythmogenesis of oscillatory patterns emerging in network motifs composed of 
inhibitory coupled tonic spiking neurons represented by the Plant model of R15 nerve cells. Such 
motifs are argued to be used as building blocks for a larger central pattern generator network 
controlling swim locomotion of sea slug Melibe leonina. 
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1. Introduction 

A plethora of vital rhythmic motor behaviors, such as heartbeat, respiratory fuuctious aud locomotiou are 
produced aud goverued by ueural uetworks called ceutral patteru geuerators (CPG s) [Selverstm , 1985; Bal 


et a/., 1988; Harder & Calabrese, 1996; Frost & Katz, 1996; Kristau et a/. [ [200^ [Katz fc Hooper , 2007| 


A CPC is a microcircuit of iuterueurous whose mutually syuergetic iuteractious autouomously geuerate 
au array of multi-phase burstiug rhythms uuderlyiug motor behaviors. There is a growiug couseusus iu 
the commuuity of ueurophysiologists aud computatioual researchers that some basic structural aud fuuc- 
tioual elemeuts must be shared by CPGs iu iuvertebrate aud vertebrate auimals. As such, we should hrst 
uuderstaud these elemeuts, hud the uuiversal priuciples, aud develop efhcieut mathematical aud computa¬ 
tioual tools for plausible aud pheuomeuological models of CPG uetworks. Pairiug experimeutal studies aud 
modeliug studies has proveu to be key to uulockiug iusights iuto operatioual aud dyuamical priuciples of 


CPGs jCilluer & Waheu, 1985; Kopell & Ermeutrout, 2004; Matsuoka[ 1987| Kopell, 1988; Cauavier et a/.. 

1994; Skiuuer et a/., 1994a; Dror et a/., 1999; Priuz et al., 2003a|. Although v 

arious circuits and models of 
IPG dynamics so robust and 

specihc CPGs have beeu c 

eveloped, it still remaius nuclear what makes the ( 

hexible [Best et a/.. 

2005; 

Belykh &; Shilnikov, 2008; Sherwood et a/., 2010; 

Koch et a/., 2011; Calabrese 

et a/., 2011; Harder 

, 2012 

. It is also unclear what mechanisms a multi-functional motor system can use 


to geuerate polyrhythmic outcomes to goveru several behaviors [Kristau , 2008; Briggmau &; Kristau, 2008 


Wojcik et a/., 2014 . Our goal is to gaiu iusight iuto the fuudameutal aud uuiversal rules goveruiug patteru 


formatiou iu complex uetworks of ueurous. To achieve this goal, we should ideutify the rules uuderlyiug 
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Fig. 1. (a) Melibe leonina lateral swim style, (b) Network bursting in swim interneurons (Si) of the Melibe swim CPG halts 

when Si3R is hyperpolarized, thus its counterpart Si3L begins tonic spiking; the photographs and in-vitro recording provided 
courtesy of A. Sakurai [Sakurai et al. 2014] 


the emergence of cooperative rhythms in simple CPG networks. 

Recently, a great deal of computational studies have been focused on a range of 3-cell motifs of bursting 
neurons coupled by chemical (inhibitory and excitatory) and electrical synapses to disclose the role of 
coupling in generating sets of coexisting rhythmic outcomes, see [Shilnikov et al , 2008; Wojcik et a/., 2011 


2014; Schwabedal et a/., 2014[ Cohens et a/.[[20T^ and references therein. These network structures 


reflect the known physiological details of various CPG networks in real animals. Next, we would like to 
explore dynamics and stability of some identifled CPG circuits constituted by 4-cells jJalil et |2013| . 
Examples of such sub-networks can be found in the cerebral crustacean stomatogastric ganglion (STG) 
jSelverstcm , 1985[ [^inz et a/.' , 2003b, 2004; Harder, 2012| , as well as in the swim CPGs of the sea slugs 


Melibe leonina (depicted during swimming in Fig. [^a)) and Dendronotus iris [Newcomb et al\ 2012 


Sakurai et al, |2011 ; [Sakurai fc Katz , 2011] . Our greater goal is to create dynamical foundations for the 


onset, morphogenesis and structural robustness of rhythmic activity patterns produced by swim CPGs in 
these animals. A pilot mathematical model of the Melibe swim CPG will be discussed in this paper. The 
circuitry shown in Fig. [^a) depicts only some core elements identifled in the biological CPG; its detailed 
diagram can be found in [Sakurai et al , 2014] . 
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Fig. 2. (a) A core circuitry of the biological Melibe swim CPG with inhibitory (•), excitatory (▼) and electrical (\/\/) synapses 

Sakurai et al. [MS] , D In-vitro voltage activity recordings from identified swim interneurons, Si2L and Si3L/R, of the Melibe 


swim CPG with the characteristic f-phase lag between the HC02 and HC03; intracellular recording provided courtesy of 


A. Sakurai [Sakurai et a/.[[T014 


Being inspired by experimental studies of voltage activity recorded from the swim CPGs of the sea 
slugs Melibe leonina and Dendronotus iris^ we would like to develop an assembly line for CPG constructors 
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Fig. 3. (a) Parabolic distribution of spike frequency within bursts produced by networked interneurons in the Melibe swim 

CPG. Recording provided courtesy of A. Sakurai and time series analysis by A. Kelley. 


made of coupled biophysically plausible models. Our first simplifying assumption is that CPGs are made 
of universal building blocks - half center oscillations (HCOs) [Hill et oL , 2003 . Loosely speaking, a HCO 
is treated as a pair of interneurons interacting with each other through reciprocally inhibitory synapses 
and exhibiting anti-phase bursting. The interneurons of a HCO can be endogenous bursters, tonic spiking 
or quiescent ones, which exhibit alternating bursting only when they inhibit each other. Theoretical stud¬ 
ies [Wang fc Rinzel , 1985 have indicated that formation of an anti-phase bursting rhythm is always based 
on slow subsystem dynamics. There are three basic mechanisms to generate alternating bursting in the 
HCO: release, escape, and post-inhibitory rebound (PIR). The first mechanism is typical for endogenously 
bursting neurons jJalil et 2010, 2012j . The other two mechanisms underlie network bursting in HCOs 
comprised of neurons, which are hyperpolarized quiescent in isolation jPerkel fc Mullon^ , 1974; Skinner 


et a/., 1994b; Angstadt et a/., 2005; Kopell & Ermentrout, 2002 . Our second assumption is that the swim 


CPG interneurons are intrinsic tonic spikers that become network bursters only when externally driven 
or coupled by inhibitory synapses, as recent experimental studies suggest [Sakurai et a/.[[^14| . The third 
assumption is that network bursting in the Melibe swim CPG is parabolie^ i.e. the spike frequency within a 
burst increases at the middle, and decreases at the ends, as one can observe from Fig. This observation 
indicates the type of neuronal models to be employed to describe network cores. Our model of choice for 
parabolic bursting is the Plant model [Plant fc Kim , 1975, 1976; Plant, 1981] . The Plant model has been 
developed to accurately describe the voltage dynamics of the R15 neuron in a mollusk Aplysia Californiea^ 
which has turned out to be an endogenous burster [Levitan fc Levitt , 1988[ . Most dynamical properties 
of the R15 neuron have been modeled and studied in detai l [Canavier et al , 1991; Bertran, 1993; Butera 
et al, 1995; Butera, 1998; Sieling & Butera, 2011; Ji et al, 2013 


2. Methods: the Plant model of parabolic bursting 

The conductance based Plant model [Plant , |1981| for the R15 neuron [Sieling fc Buteraj |2011[ located in 
the abdominal ganglion of a slug Aplysia Californiea is given by the following set of ordinary differential 
equations derived within the framework of the Hodgkin-Huxley formalism to describe the dynamics of the 
fast inward sodium [Na], outward potassium [K], slow TTX-resistant calcium [Ca] and an outward calcium 
sensitive potassium [KCa] currents: 


CmV ~ ^Na ICa I KCa I leak I ext 


syu' 


( 1 ) 


The last three currents are the generic ohmic leak Iieak^ external constant I^xt and synaptic Isyn currents 
flowing from a pre-synaptic neuron. The full details of the representation of the currents employed in the 
model are given in the Appendix below. 
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Fig. 4. (a) Endogenous bursting in the Plant model as alternations of tonic spiking activity and quiescent periods, (b) Single 

burst featuring a characteristic spike frequency increase in the middle of each burst, (c) Parabolic shape of the frequency 
distribution of spikes within a burst is a feature of this kind of bursting. The parameters are p = 0.00015ms“^, Kc = 0.00425ms“^ 
and Tx = 9400?ns. 


There are two bifurcation parameters in the individual model. The first one is the constant external 
current, Iext-> which is set I^xt - 0. Following [Shilnikov , 2012] , the other bifurcation parameter, A, is 
introduced in the slowest equation: 

Ca = p {K.xiVca -V^/\)-Ca) (2) 

describing the concentration of the intracellular calcium in the Plant model. By construction, A is a 
deviation from a mean value of the reversal potential Vca - IdOmF evaluated experimentally for the 
calcium current in the R15 cells. As such, this makes A a bifurcation parameter. Secondly its variations 
are not supposed to alter the topology of the slow motion manifolds in the 5D phase space, which are called 
tonic spiking and quiescent in the mathematical neuroscience context, as they are made of, respectively, 
round periodic orbits and equilibrium states [of the slow subsystem] of the mod el (Fig. [^. 

At A = 0, the neuron is an endogenous burster, see Fig. According to [Rinzel fc Lee , |1987j , this 
type of bursting is termed parabolic. The reason for this term is that the spike frequency within bursts is 
maximized in the middle of bursts and minimized at the beginning and the end (see Fig. [^). The parabolic 
structure of a burst is due to the calcium-activated potassium current. Its magnitude is determined by 
the intracellular calcium concentration. As the intracellular calcium concentration increases, the calcium 
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Fig. 5. Bursting (green) orbit recursively switching between two slow-motion critical manifolds: tonic spiking, with 

a characteristic fold and originating through a sub-critical Andronov-Hopf (AH) bifurcation from a depolarized equilibrium 
state, and quiescent, Meq (orange curve), projected onto the (h,P) and slow Ca variables of the of the Plant model; a plane 
represents the synaptic threshold, ©s^/n = OmP. 
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Fig. 6. Responses of the bursting neuron (A = OmV) on the synaptic drive Isyn = gsyn{V - Vrev)- (a) Excitatory synaptic 
drive with gsyn = 0.002n*S' and Vrev = ^OmV applied at t = 80sec switches the neuron from bursting to tonic spiking activity, (b) 
The inhibitory drive with gsyn = 0.005nS and Vsyn = -80mV halts bursting and makes the neuron hyperpolarized quiescent. 


dependent potassium current gets activated, which causes an increase of the inward potassium current. As 
the membrane potential increases over a threshold value, the intracellular calcium concentration decreases, 
as well as the inward potassium current (see Eq. (20) in the Appendix). The parabolic distribution of 
spikes within bursts is shown in Fig. The instant frequency value is calculated by the reciprocal of each 
inter-spike interval. Panels b and c of Fig. [^clearly disclose the parabolic inter-spike structure of bursts. 

It was shown in [Rinzel fc Lee , [T987| that the mechanism underlying a transition between quiescent and 
tonic spiking of bursting in the Plant model is due to a homoclinic bifurcation of a saddle-node equilibrium 
state [Shilnikov , 1963; Afraimovich et a/., 2014| . This bifurcation occurs in the fast 3D (F,/i,n)-subspace 
of the model and is modulated by the 2D slow dynamics in the ((7a, x)-variables, which are determined 
by slow oscillations of the intracellular calcium concentration [Plant fc Kim , 1975, 1976j . The unfolding 
of this codimension-one bifurcation includes an onset of a stable equilibrium, which is associated with a 
hyperpolarized phase of bursting, and on the other end, an emergent stable periodic orbit that is associated 
with tonic spiking phase of bursting. The period of this stable orbit decreases, as it moves further away 
from the saddle-node equilibrium mediated by decreasing calcium concentration. The period of the tonic 
spiking orbit grows with no upper bound as it approaches the homoclinic loop of the saddle-node [Shilnikov 


et a/., 1998,2001 


Variations of A change the duty cycle of bursting, which is a ratio of the active tonic spiking phase of 
bursting to its period. Decreasing A reduces the inactive, quiescent phase of bursting, i.e. increases its duty 
cycle. Zero duty cycle is associated with the homoclinic saddle-node bifurcation that makes the neuron 
hyperpolarized quiescent. This corresponds to an emergence of stable equilibrium state for all dynamical 
variables of the model Q. In other words, decreasing A makes the active phase longer, so that below a 
threshold A = -32mV the neuron switches to tonic spiking activity. Tonic spiking activity is associated with 
the emergence of a stable periodic orbit in the fast (F, /i, n)-subspace, while the ((7a, x)-variables of the slow 
subspace converge to a stable equilibrium state. As such, bursting occurs in the Plant and similar models 
due to relaxation of periodic oscillations in the 2D ((7a, x)-subspace, which slowly modulates fast tonic 
spiking oscillations in the (F, /i,n) variables. The relaxation limit cycles emerge from one and collapse into 
the other equilibrium state in the ((7a, x)-plane through Andronov-Hopf bifurcations, which can be sub- 
or super-critical. At the transitions between bursting and tonic spiking, and bursting and hyperpolarized 
quiescence, the neuron can produce chaotic dynamics, which are basically due to the membrane potential 
oscillatory perturbations of plain canards at the folds of the relaxation cycle. 
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Fig. 7. Tonic spiking neuron 1 at A = -34mV near the bifurcation transition between tonic spiking and bursting is forced to 
become a network burster with an application of an inhibitory drive with gl^n = O.OOlnS, from the pre-synaptic neuron 2 at 
t = GOsec. Halting the inhibitory drive restores tonic spiking activity in the targeted neuron (not shown). 


3. Endogenous and network bursting. Inhibitory and excitatory drives 

A half-center oscillator is a network of two neurons coupled by reciprocally inhibitory synapses that robustly 
produces bursting in alternation, or anti-phase bursting. Such a network can be multistable, i.e. produce 
other bursting rhythms as well, such as synchronous bursting [Jalil et 2010| and r hythmic outcomes 
with slightly shifted phase lags between the endogenously bursting neurons [Jalil et 2012 . 

In this study, the synaptic current Igyn is modeled through the fast threshold modulation (FTM) 


approach [Kopell fc Somers' , 1993j . The synapses are assumed to be fast and non-delayed, which is true for 
the swim CPG in both sea slugs under consideration. The synaptic current is given by 


1 


Isyn - gsyn{Vpost Esyn) ^ ^ ^-kiVpre-Qsyn) 


( 3 ) 


where gsyn is the maximal conductance of the current, which is used as a bifurcation parameter of the 
networked model; Vpost{t) and Vpre{t) are the voltages on the post-synaptic (driven) and pre-synaptic 
(driving) neurons; Egy^ is the synaptic reversal potential. To make Isyn excitatory, we set Egyn = dOmF, 
while in the inhibitory case we set Egyn = -SOml/. In Eq. ([^, the second term is a Boltzmann coupling 
function that quickly, {k - 100), turns the synaptic current on and off as soon the voltage, Vpre^ of the 
(driving) pre-synaptic cell(s) raises above and falls below the synaptic threshold, here Qgyn - OmE (Fig.j^. 


To model the constant synaptic drive onto the post-synaptic neuron, we assume that Vpre > © 


^syn- 


This allows us to calibrate the state of the post-synaptic neuron, and to determine the drive threshold 
that separates the qualitatively distinct states of the individual and networked neurons. This statement is 
illustrated in Fig. by simulating responses of the endogenous parabolic burster to network perturbation. 
Figure a) shows, with a properly adjusted excitatory drive, that the endogenous burster switches into 
tonic spiking activity. On the other hand, bursting in the networked neuron can be halted when it receives 
a sufficient inhibitory drive from the pre-synaptic neuron of the network (Figure [^b)). Eliminating either 
drive makes the post-synaptic neuron return to its natural state, i.e. these experiments de-facto prove that 
the neuron is mono-stable for the given parameter values. 

A HCO, in the canonical Brown definition [Graham-Brown , 1911[ , is a pair of neurons bursting in 
anti-phase when they are networked by inhibitory synapses. In isolation, such neurons are not endogenous 
bursters but tonic spikers instead, or remain quiescent [Harder fc Calabrese , 1996] . There are multiple 
mechanisms underlying such anti-phase bursting, or, more accurately, anti-phase oscillations in HCOs and 
CPGs made of relaxation oscillators [Kopell fc Ermentront , 2002; Daun et a/., 2009] . The list includes the 
well studied mechanisms of post-inhibitory rebound and escape for quies cent neurons [Perkel fc Mnlloney 


1974; Wang &; Rinzel, 1985; Skinner et a/., 1994b; Destexhe et a/., 1994; Matveev et a/., 2007 , as well as 


less-known mechanisms of HCOs constituted by intrinsically spiking neurons. Such networks utilizing the 
Plant models are discussed below. 

To construct such a HCO with relatively weak inhibitory coupling, the Plant model must be first set 
into the tonic spiking mode. This is done by setting the bifurcation parameter, A = -34mK, see Fig. 
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Fig. 8. Anti-phase network bursting produced by a HCO of two Plant neurons as soon as the inhibition is turned on. Blocking 
the inhibition restores tonic spiking activity in both neurons, and vice versa. Here, the the network parameters are = 
O.OOSnS and Esyn = -SOmV, and the parameters of the individual neurons are the following: A = -60mV, p = O.OOOSms”^, 
Kc = 0.0085ms“\ Tx = 235ms and XociV) = 1/(1 + 


Next, we consider a unidirectional network where the tonic spiking neuron 1 starts receiving, an inhibitory 
drive of ggyn = O-OOlnS" from the post-synaptic neuron 2 at t = GOsec. The inhibitory drive is sufficient to 
shift the post-inhibitory neuron over the bifurcation transition back into bursting activity. The minimal 
inhibitory drive must be increased proportionally to make the targeted neuron a network burster whenever 
it stays further away from the bifurcation transition between tonic spiking and bursting in isolation. 


4. Forming a half-center oscillator 

In this section, we discuss the dynamics of half-center oscillators made of two tonically spiking Plant neurons 
reciprocally coupled by inhibitory synapses. As before, we describe such synapses within the framework 
of the fast threshold modulation (FTM) paradigm using Eq. to match the shape and magnitude of 
inhibitory postsynaptic potentials (IPSPs) in the post-synaptic neurons. IPSPs are the indicators of the 
type and the strength of synapses in the network. 

We perform simulations in a fashion that is analogous to the dynamic clamp technique used in neu¬ 
rophysiological experiments. The approach involves the dynamic block, restoration and modulation of 
synaptic connections during simulation. These modeling perturbations should closely resemble the exper¬ 
imental techniques of drug-induced synaptic blockade, modulation, wash-out, etc. Restoring the chemical 
synapses during a simulation makes the HCO regain network bursting activity with specific phase charac¬ 
teristics. Depending on the coupling strength as well as the way the tonically spiking neurons are clamped, 
the network bursting may change phase-locked states, i.e. be potentially multi-stable. Experimental ob¬ 
servations also suggest specific constraints on the range of coupling strengths of the reciprocal inhibition, 
such that the networks stably and generically achieve the desired phase-locking. 

Figure [^demonstrates the stages of anti-phase bursting formation in the HCO. The uncoupled neurons 
are initiated in tonic spiking mode. After turning on the reciprocally inhibitory synapses ggyn = 0.008n5, 
the HCO quickly transitions to the regime of robust anti-phase bursting. Turning off the synapses restores 
the native tonic spiking activity in both neurons. Turning on the reciprocal synapses makes the HCO 
regain the network bursting. Note that the length of transients from tonic spiking to network bursting 
depends on the strength of the synaptic coupling for the fixed parameters of the individual Plant neurons. 
By comparing the magnitude of IPSPs in the voltage traces represented in Figs. and one can conclude 
that the coupling in the later case is weaker. This is why, the onset of network bursting in the HCO is less 
pronounced. 

Our modeling studies agree well with experimental recordings from the identified interneurons in the 
Melibe swim CPC which suggests that the observed bursting is due to synergetic interactions of interneurons 
of the network [Saknrai et a/. , 2014| . One can see from Fig. [^b) that network bursting in the biological 
HCO formed by two Si3 interneurons of the Melibe swim CPC is seized as soon as the right one, Si3R, 
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Fig. 9. Onset of emergent network anti-phase bursting in the HCO with reciprocally inhibitory. [Esyn = -80mV) synapses 
at giyn = 0.0073n>S. 


receives a negative current pulse that makes it hyperpolarized quiescent, while its left bursting counterpart, 
Si3L, turns into tonic spiking activity instead. Moreover, one can deduct from the wiring diagram of the 
CPG depicted in Fig. [^a) and the analysis of voltage traces represented in Fig. [^b) that the interneuron 
Si2L becomes a tonic-spiker as soon as the pre-synaptic interneuron Si3R stops inhibiting it (compare 
with Fig. 1^) This further supports the assertion that the swim CPG is made of intrinsically tonic spiking 
interneurons. 

To test the robustness of network anti-phase bursting to perturbation and to calibrate the necessary 
influx of reciprocal inhibition generated by the Plant neurons, we consider a HCO with excitatory autapses. 
The objective here is to determine an equivalent amount of excitatory drive to be projected onto the post- 
inhibitory network burster to cancel out the inhibitory drive and shift it back to the initial tonic-spiking 
mode. 

An autapse is a synapse of a neuron onto itself, where the axon of the neuron ends on its own dendrite. 
After their discovery [Van der Loos fc Glaser , 1972| autapses have been observed in a range of nervous 
systems. The autapses are arguably to be responsible for tuning of neural networks. This particular con¬ 
figuration of the HCO depicted in Fig. [^is formally motivated by the swim CPG circuitry, see Fig. [^a). 
One can see from it that the interneurons of the bottom HCO receive excitatory drives from the top in¬ 
terneurons forming the top HCO. We would like to find the threshold over which the neurons no longer 
form a stably bursting HCO. This would allow us to calibrate and quantify the relative strengths of the 
mixed synaptic connections in the swim CPG models. 

In this HCO configuration, each neuron inhibits its counterpart and self-excites through the autapse. 
Both autapses are introduced to the model using the FTM approach with Eaut = 40mH. In this experiment. 




Time(sec) 


Fig. 10. Turning on the excitatory autapses at = O.OlGnS in the HCO with gl^n = 0.0073nS halts pronounced network 
bursting. 
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Fig. 11. Assembly line of the Melibe swim CPG model out of four intrinsically tonic spiking Plant neurons. First the reciprocal 
inhibition between Si3R and Si3L is turned on, followed by turning on the reciprocal inhibition between Si2R and Si2L, and next 
simultaneous turning on unidirectional cross-lateral inhibition from Si3R(L) projected onto Si2L(R), and bi-lateral excitation 
originating from Si2R(L) down onto Si3R(L). After a short transient, the CPG model exhibits the desired 3/4 phase shift lag 
between Si2L and Si3L. Compare with voltage traces of the biological CPG in Fig.[^b). 


the conductance values for inhibitory synapses are set at = 0.0073n5. This is sufficient for the HCO to 
generate robust anti-phase bursting as seen in Fig.[^ Next, we add the autapses along with inhibition and 
gradually increase g^. We found that increasing g^ proportionally increase the delay. At g^ = 0.016n5, 
the network stops exhibiting anti-phase bursting. We note that unlike a permanent excitatory drive from 
pre-synaptic neurons, an introduction of the excitatory autapse, acting only when the self-driving neuron 
is above the synaptic threshold, is effectively perturbation equivalent for the calibration purpose. 


5. Assembly line of a Melibe swim CPG 

In this final section, we put together a pilot model of the Melibe swim CPG according to a circuitry based 
on identified interneurons and synapses; its wiring diagram is sketched in Fig. [^a). This network model 
is made of the two HCOs constituted by tonic spiking Plant neurons. We would like to find out whether 
this sample CPG model can already produce phase lags similar to those between bursting interneurons in 
the biological CPG. For the sake of simplicity, we do not include Si4R/L interneurons in the model and 
we also omit electrical synapses. It is known from experimental studies [Saknrai et al , 2014 that blocking 
chemical, inhibitory and excitatory synapses between the interneurons may be sufficient to break down 
the motor pattern by the network. Figure [^b) points out that the interneurons of either HCO burst in 
anti-phase and there is the characteristic 3/4 phase lag between the burst initiation in the neurons Si2L 
and Si3L, as well as between Si2R and Si3L. This phase lag is repeatedly observed in both adult and 
juvenile animals. 

As before, we use the Plant neurons initiated in the tonic spiking mode, relatively close to the transition 
to bursting. Initial conditions of the neurons are randomized. After letting the neurons settled down to tonic 
spiking activity, the network connections are turned on. As Fig. shows, with the reciprocal inhibition 
being first turned on, the bottom interneurons Si3L and Si3R become anti-phase network bursters, and 
so do Si2R and Si2L as soon as the reciprocal inhibition between them is turned them on, too. At this 
stage, the CPG model is formed by two uncoupled HCOs. A few seconds later, they become coupled by 




















































































































































10 Alagam and Shilnikov 


simultaneous turning on the unidirectional cross-lateral inhibition from Si3R(L) projected onto Si2L(R), 
and bi-lateral excitation from Si2R(L) down onto Si3R(L). One can see from this figure that all four 
interneurons of the CPG model exhibit network bursting with the desired phase lags. These are 0.5 (half 
period) between the interneurons of each HCO, and 3/4 (a fraction of the network period) between the 
HCOs, or between the corresponding reference interneurons Si2L and Si3L. We note that such a phase 
shift was reported in a similar Melibe swim CPG constituted by endogenous bursters; that model also 
incorporated electrical synapses [Jalil et oT , 2013| . There is a great room for improvement of CPG network 
models to include other identified interneurons and to incorporate additional electrical synapses to find 
out whether additions of new elements can stabilize or desynchronize the desired bursting pattern as it 
was done using the Poincare return maps for endogenous bursters [Wojcik et a/. , 2014 . Of our special 
interest are various problems concerning structural stability of the network, and its robustness (Lyapunov 
stability) for bursting outcomes subjected to perturbations by pulses of the external current, as well as 
reductions to return maps between burst initiations in constituent neurons. These questions are beyond 
the scopes of the given examination and will be addressed in full detail in our forthcoming publications 
soon. The question about a possible linking of the characteristic | phase lag and the Melibe leonina lateral 
swim style is the paramount one among them. 


6. Summary 

We have discussed a basic procedure for building network bursting CPGs made of intrinsically tonic spiking 
neurons. As a model for such networks, we have employed the biophysically plausible Plant model that 
was originally proposed to describe endogenous bursting R15-cells in the Aplysia mollusk. Such bursting 
was intracellularly recorded, and identified as parabolic, from the known interneurons in the swim CPGs of 
two sea slugs: Melibe leonina and Dendronotus iris. There is experimental evidence that bursting in these 
swim CPGs is due to synergetic interactions of all constituent neurons that are intrinsic tonic-spikers in 
isolation. To model the Melibe swim CPG, we have first examined dynamical and structural properties of 
the Plant model and its responses to perturbations. These perturbations include inhibitory and excitatory 
inputs from pre-synaptic neurons in the network. We have identified the transition boundary beyond 
which the bursting Plant model becomes a tonic-spiker and shifted it slightly over the threshold using an 
introduced bifurcation parameter. We have shown that the perturbed/calibrated Plant neuron, exhibiting 
intrinsically tonic spiking activity, becomes a network burster when it receives an inhibitory drive from 
a pre-synaptic neuron. By combining two such neurons, we have created a genuine half-center oscillator 
robustly producing anti-phase bursting dynamics. We have also considered a HCO configuration with two 
excitatory autapses to assess the robustness of anti-phase bursting with respect to excitatory perturbations. 
Finally, we have employed all necessary components to assemble a truncated model of the Melibe swim CPG 
with the characteristic 3/4-phase lags between the bursting onsets in the four constituent interneurons. 
In future studies, we plan to examine the dynamics of the CPG models with all synaptic connections, 
including electrical, as well as incorporating additional identified interneurons. We will also explore their 
structural stability, robustness and potential multi-stability of their bursting outcomes with various phase 
lags. An additional goal is to find out whether the motor pattern with the 3/4-phase lags will persist 
in networks with interneurons represented by other mathematical models including phenomenologically 
reduced ones. Potentially, these findings shall provide a systematic basis for comprehension of plausible 
biophysical mechanisms for the origination and regulation of rhythmic patterns generated by various CPGs. 
Our goal is to extend and generalize the dynamical principles disclosed in the considered networks for other 
neural systems besides locomotion, such as olfactory cellular networks. 
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8. Appendix: the conductance based Plant model 


The model in this study is adopted from [Plant 
governed by the following equation: 


1981j . The dynamics of the 


membrane potential, is 


CmV = -iNa - Ik - ICa “ IxCa “ heak “ Isyn, (4) 

where Cm = is the membrane capacitance, I^a is the Na'^ current, Ik is the current, Ica 

is the Ca'^‘^ current, IxCa is the Co?^ activated current, Iieak is the leak current, Igyn is the synaptic 
current. The fast inward sodium current is given by 


iNa = “ Vno), 


( 5 ) 


where the reversal potential V^a - SOmF and the maximum Na^ conductance value gNa - 4n*S. The 
instantaneous activation variable is defined as 



m (V) - “w(^) 

^ am{V)+P^{Vy 

(6) 

where 

“ hxp((50 - m}l0) - 1 • 

( 7 ) 

while the dynamics of inactivation variable h is given by 



• hoo{V)-h 

ThiV) ’ 

(8) 

where 

hoo{V)= X and ThiV) = , , ,, 

( 9 ) 

with 

MV) = 0.07exp((25-r,)/20) apd «»") %^p ((55 _ ’ 

(10) 

where 

,, 1271/+ 8265 ,, 

14 =- mV. 

105 

(11) 

The fast potassium current is given by the equation 



lK = 9Kn\V-VK), 

(12) 


where the reversal potential is Vk = -I^mV and the maximum conductance value is qk = 0.3n*S.The 
dynamics of inactivation gating variable is described by 


nooiV)-n 

’ 

where 

^ ~ ah{v) + ph{v) """" ^ ~ c^hiv)+Hhivy 

with 

55 - y 

aniV) = 0m - y and /3„(y) = 0.125exp((45 - V)/80). 

exp((55 - Vs)/10) - 1 

The TTX-resistant calcium current is given by 

ICa = 9CaX{V - Vca), 


(13) 


(14) 


(15) 

(16) 





















REFERENCES 13 


where the reversal potential is Vca - 140my and the maximum Co?^ conductance is gca - 0.03n*S. The 
dynamics of the slow activation variable is described by 


Xoo(V)-X 

X = 


where 


Xoo{V) 


exp(-0.3(1/ + 40)) + l 


r.iV) ’ 

and Tx{V) = 9400m5. 


The outward activated current is given by 


IxCa = RKCa ^ } , i 


(17) 

(18) 

(19) 


0.5 + [Ca]^ 

where the reversal potential is Vca = l^OmV. The dynamics of intracellular calcium concentration is 
governed by 

Ca = p [K,x{Vca -V)- [Ca],], (20) 

where the reversal potential is Vca = 140my, and the constant values are p = 0.00015my^ and Kc = 
0.00425my“^. The leak current is given by 

Ileak=gL{V-VL)^ ( 21 ) 

where the reversal potential Vl = -40mV and the maximum conductance value = 0.000377.5. The synaptic 
current is defined as 


9syni^post ^rev^ 


( 22 ) 


2 _|_ ^—k(^Vpre~^syn) 

with the synaptic reversal potential Vpost = -80mV for inhibitory synapses and Vpost = ^OmV for excitatory 
synapses and the synaptic threshold Qsyn = OrnV^ and k = 100. 
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